--------------------------------------------------------------------------------
-- This package is a thick Ada binding to the Multiple Precision Floating Point
-- Reliable Library. See <http://www.mpfr.org/>.
--
-- Copyright (C) 2009 Vincent DIEMUNSCH.
-- GNU General Public License.
--
-- This program is free software: you can redistribute it and/or modify it under
-- the terms of the GNU General Public License as published by the Free Software
-- Foundation, either version 3 of the License, or (at your option) any later
-- version.
--
-- This file is distributed in the hope that it will be useful, but WITHOUT ANY
-- WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
-- A PARTICULAR PURPOSE. See the GNU General Public License for more details.
--
-- You should have received a copy of the GNU General Public License along with
-- this program.  If not, see <http://www.gnu.org/licenses/>.
--------------------------------------------------------------------------------
--
-- This package declares a new private type : MPFR_Float, that is an abstraction
-- of a floating point value of arbitrary precision.
--
-- This packages uses a new approach to dynamically adjust the precision of a
-- MPFR_Float, everytime a value is computed (and not, as mpfr does, at the
-- initialization of the mpfr_t variable). The goal is to try to have a precision
-- that is of the same order of magnitude as the accuracy.
-- Before any operation, the precision required for the result is set according to
-- the following rules:
--    1) if the result is exact, then precision is the default precision
--       for instance : , "1/3", "pi/4", "10.5", "1.45e-3".
--    2) for a unary operation, the precision of the result is the precision
--       of the operand.
--       for instance : "abs(x)", "1/x", "x+n" where n is an integer.
--    3) for a binary multipliying operation the precision of the result is the
--       smallest precision of the two operands.
--       for instance : "x*y", "x mod y", etc.
--    4) for a binary adding operation, the precision of the result is the
--       smallest precision of :
--        * the greatest operand in absolute value,
--        * the smallest operand in absolute value, adjusted to the difference
--          of exponants.
--       for instance : "1.111111 + 1.0e-5" (resp. a precision of 7 digits and
--                       2 digits) will give "1.111121" (7 digits = 2 - (0-5))
--                       and not "1.1" (2 digits).
--
-- If, an adjustement of precision is required, in a computation, for instance
-- when the result is computed through the use of a convergent sequence, one
-- can use the functions "Unchecked_Pecision_Adjustement", with 2 variants :
--   1) an absolute precision expressed in digits
--   2) a relative precision, throug the comparison to another value :
--      new precision = exponant(x) - exponant(x - value) + n digits (n > 3).
--      for instance :
--          u(n) = 1/u(n-1) + 1   (Computing the golden number)
--          u(n) = 1.615 (4 digits) and u(n-1) = 1.625, error : 1e-2,
--          then the precision will be set to 0 - (-2) + n digits, at least
--          5 digits ("1.6150")


with Ada.Finalization;
with GMP.Integers;

package MPFR.Floats is

   type MPFR_Float is private;

   -- Private type so that all the following functions are NOT primitive
   -- operations.

   -- PRECISION : --------------------------------------------------------------

   procedure Set_Default_Precision (Decimal_Digits : Positive);
   -- Set the default precision whenever it is used, according to ARM G.2.2.
   -- Note : if you want to precisely set this precision in bits simply use the
   --        procedure mpfr_set_default_prec (prec : mpfr_prec_t) from the
   --        MPFR package.
   function Get_Default_Precision return Positive;

   procedure Unchecked_Precision_Adjustement
     (x              : in out MPFR_Float;
      Decimal_Digits : Positive);

   procedure Unchecked_Precision_Adjustement
     (x              : in out MPFR_Float;
      Error          : in     MPFR_Float;
      Decimal_Digits : Positive := 3);

   -- Return the value of x, either rounded or extended with trailing zeros.
   -- Adjust the precision either to the absolute precision given in digits, or
   -- by adjusting it to the given digits plus the order of magnitude of the
   -- the difference between x and the reference value.

   type Rounding_Mode_Kind is (Round_to_Nearest,
                               Round_Toward_Zero,
                               Round_Toward_Plus_Infinity,
                               Round_Toward_Minus_Infinity);
   -- we declare here exactly the same type as the mpfr_rnd_t from the
   -- MPFR package, but for the sake of clarity and simplicity of the user.

   procedure Set_Default_Rounding_Mode (Rounding : Rounding_Mode_Kind);
   -- Set the default rounding used by this package, every time an operation
   -- is done on a MPFR_Float. The default rounding mode is
   -- to nearest initially.

   procedure Set_Precision_Check (Decimal_Digits : Positive);
   procedure Inhibate_Precision_Check;
   -- These functions enable or inhibate the checking during assignement to a
   -- variable, that the value has a precision above the curret minimal required
   -- precision.
   -- If needed, the Not_Accurate exception is raised.

   Default_Precision : constant Natural := 0;



   -- INITIALISATION : ---------------------------------------------------------
   function To_MPFR_Float (Image          : String;
                           Decimal_Digits : Natural := Default_Precision)
                           return MPFR_Float;

   -- This function creates Floating Point from a value given by a string
   -- that must conform to the Ada Standard for Float literals (decimal literals
   -- and base literals), and whose precision is at least the number of
   -- decimal_digits given, or if this value is nul the Default_Precision.
   -- The base can be above 16, but it must be under (or equal to) 36 (and in
   -- that case we use digits from 0 .. 9, A .. Z).
   -- We tolerate a negative value (i.e. the '-' at the beginning), even if in
   -- Ada the '-' is seen as a unary operator and not part of the literal.
   -- Exemples :
   --      Decimal literals : 12.0   0.0   0.456  3.14159_26
   --      Based   literals : 16#F.FF#E+2   2#1.1111_1111_1110#E11  = 4095.0

   Data_Error : Exception;

   subtype Base_Range is Integer range 2 .. 36;

   function Image (Item           : MPFR_Float;
                   Base           : Base_Range := 10;
				   Full_Mantissa  : Boolean :=  False)   return String;
   -- This function return a Float literal (decimal literal or base literal)
   -- that conform to the Ada Standard except the '-' at the beginning and the
   -- fact that base can be any integer between 2 and 36.
   -- Note that if Full_Mantissa is true, mfr will print the full content of the
   -- Mantissa and not only the part required by the decimal precision according
   -- to ARM G.2.2, S'Model_Mantissa


   function Not_A_Number   return MPFR_Float;
   function Plus_Infinity  return MPFR_Float;
   function Minus_Infinity return MPFR_Float;

   function Pi return MPFR_Float;
   -- Pi at the current default precision.

   ----------------------------
   -- "Standard" description --
   ----------------------------

   -- (Mimics Standard package)
   -- The predefined operators for this type are as follows:

   function "="   (Left, Right : MPFR_Float) return Boolean;
   function "<"   (Left, Right : MPFR_Float) return Boolean;
   function "<="  (Left, Right : MPFR_Float) return Boolean;
   function ">"   (Left, Right : MPFR_Float) return Boolean;
   function ">="  (Left, Right : MPFR_Float) return Boolean;
   pragma Inline ("=", "<", "<=", ">", ">=");


   function "abs" (Right : MPFR_Float) return MPFR_Float;
   function "-"   (Right : MPFR_Float) return MPFR_Float;
   pragma Inline ("abs", "-");

   function "+"   (Left, Right : MPFR_Float) return MPFR_Float;
   function "-"   (Left, Right : MPFR_Float) return MPFR_Float;
   function "*"   (Left, Right : MPFR_Float) return MPFR_Float;
   function "/"   (Left, Right : MPFR_Float) return MPFR_Float;

   function "**"  (Left : MPFR_Float; Right : Integer'Base) return MPFR_Float;


   ----------------------------
   -- With Standard Integers --
   ----------------------------

   function To_MPFR_Float (Value          : Integer'Base;
                           Decimal_Digits : Natural := Default_Precision)
                           return MPFR_Float;

   -- This function creates a Floating Point value from the given integer value
   -- whose precision is at least the number of decimal_digits given, or if this
   -- value is nul the Default_Precision.


   -- In addition, the following operators are predefined for the root
   -- numeric types:
   function "*" (Left : MPFR_Float;   Right : Integer'Base) return MPFR_Float;
   function "*" (Left : Integer'Base; Right : MPFR_Float)   return MPFR_Float;
   function "/" (Left : MPFR_Float;   Right : Integer'Base) return MPFR_Float;

   -- and we add :
   function "+" (Left : MPFR_Float;   Right : Integer'Base) return MPFR_Float;
   function "+" (Left : Integer'Base; Right : MPFR_Float)   return MPFR_Float;
   function "-" (Left : MPFR_Float;   Right : Integer'Base) return MPFR_Float;
   function "-" (Left : Integer'Base; Right : MPFR_Float)   return MPFR_Float;
   function "/" (Left : Integer'Base; Right : MPFR_Float)   return MPFR_Float;

   --------------------------
   -- with Standard Floats --
   --------------------------
   generic
      type F is digits <>;
   function Generic_To_MPFR_Float (Value : F) return MPFR_Float;

   function To_MPFR_Float (Value          : Float;
                           Decimal_Digits : Natural := Default_Precision)
                           return MPFR_Float;

   function To_MPFR_Float (Value          : Long_Float;
                           Decimal_Digits : Natural := Default_Precision)
                           return MPFR_Float;
   pragma Inline (To_MPFR_Float);


   generic
      type F is digits <>;
   function Generic_Round (X : MPFR_Float) return F;

   function Round (X : MPFR_Float) return Float;
   function Round (X : MPFR_Float) return Long_Float;
   function Round (X : MPFR_Float) return Long_Long_Float;


   -- This functions compares the multiprecision value with the [long_]float value
   -- and returns the error value = Value - Exact_Value.
   function Compare (Value : Float;      Exact_Value : MPFR_Float) return Float;
   function Compare (Value : Long_Float; Exact_Value : MPFR_Float) return Long_Float;


   -----------------------------
   -- with Unbounded_Integers --
   -----------------------------

   subtype Unbounded_Integer is GMP.Integers.Unbounded_Integer;

   Best_Precision : constant Natural := 1;

   function To_MPFR_Float (Value          : Unbounded_Integer;
                           Decimal_Digits : Natural := Default_Precision)
                           return MPFR_Float;
   -- This function creates A Floating Point value from the given integer value
   -- whose precision is at least the number of Decimal_Digits given, or
   --  * if this value is nul the Default_Precision,
   --  * if this value is   1 the max of the Default Precision and the number
   --    of decimal digits of the unbounded integer

   function Round_To_Unbounded_Integer(Item     : MPFR_Float)
                                       return Unbounded_Integer;
   function Round_To_Unbounded_Integer(Item     : MPFR_Float;
                                       Rounding : Rounding_Mode_Kind)
                                       return Unbounded_Integer;
   -- Round to an unbounded integer, using either the default rounding mode
   -- or the specified mode.


   --------------------------
   -- Elementary functions --
   --------------------------

   function Sqrt (X    : MPFR_Float)        return MPFR_Float;
   function Log  (X    : MPFR_Float)        return MPFR_Float;
   function Log  (X    : MPFR_Float;
                  Base : Integer'Base)      return MPFR_Float;
   function Exp  (X    : MPFR_Float)        return MPFR_Float;
   function "**" (Left : MPFR_Float;
                  Right: Unbounded_Integer) return MPFR_Float;
   function "**" (Left, Right : MPFR_Float) return MPFR_Float;

   function Sin    (X : MPFR_Float) return MPFR_Float;
   function Cos    (X : MPFR_Float) return MPFR_Float;
   function Tan    (X : MPFR_Float) return MPFR_Float;
   function Sec    (X : MPFR_Float) return MPFR_Float;
   function Csc    (X : MPFR_Float) return MPFR_Float;
   function Cot    (X : MPFR_Float) return MPFR_Float;

   function Arcsin (X : MPFR_Float) return MPFR_Float;
   function Arccos (X : MPFR_Float) return MPFR_Float;
   function Arctan (X : MPFR_Float) return MPFR_Float;

   function Sinh   (X : MPFR_Float) return MPFR_Float;
   function Cosh   (X : MPFR_Float) return MPFR_Float;
   function Tanh   (X : MPFR_Float) return MPFR_Float;
   function Sech   (X : MPFR_Float) return MPFR_Float;
   function Csch   (X : MPFR_Float) return MPFR_Float;
   function Coth   (X : MPFR_Float) return MPFR_Float;
   function Arcsinh(X : MPFR_Float) return MPFR_Float;
   function Arccosh(X : MPFR_Float) return MPFR_Float;
   function Arctanh(X : MPFR_Float) return MPFR_Float;

   function Factorial (N : Natural'Base) return MPFR_Float;
   function Gamma     (X : MPFR_Float)   return MPFR_Float;

   ------------------
   -- Exceptions : --
   ------------------

   MPFR_Underflow    : exception;
   -- occurs when the exact result of a function is a non-zero real number and
   -- the result obtained after the rounding, assuming an unbounded exponent
   -- range (for the rounding), has an exponent smaller than the minimum exponent
   -- of the current range. In the round-to-nearest mode, the halfway case is
   -- rounded toward zero.
   MPFR_Overflow     : exception;
   -- occurs when the exact result of a function is a non-zero real number and
   -- the result obtained after the rounding, assuming an unbounded exponent
   -- range (for the rounding), has an exponent larger than the maximum exponent
   -- of the current range. In the round-to-nearest mode, the result is infinite.
   MPFR_Not_a_Number : exception;
   -- occurs when the result of a function is a NaN.
   MPFR_Inexact      : exception;
   -- occurs when the result of a function cannot be represented exactly and
   -- must be rounded.
   MPFR_Range_Error  : exception;
   -- occurs when a function that does not return a MPFR number (such as
   -- comparisons and conversions to an integer) has an invalid result.

   Not_Accurate : exception;
   -- This exception might be raised by the "Precision_Check".
   -- It is not a MPFR exception in fact !


private

   type MPFR_Float is new Ada.Finalization.Controlled with record
      Value    : MPFR.mpfr_t;
   end record;

   overriding procedure Initialize (Object: in out MPFR_Float);
   overriding procedure Adjust     (Object: in out MPFR_Float);
   overriding procedure Finalize   (Object: in out MPFR_Float);

end MPFR.Floats;

